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ABSTRACT 

We construct dynamical models of the Milky Way’s Box/Peanut (B/P) bulge, using the re¬ 
cently measured 3D density of Red Clump Giants (RCGs) as well as kinematic data from the 
BRAVA survey. We match these data using the NMAGIC Made-to-Measure method, starting 
with N-body models for barred discs in different dark matter haloes. We determine the total 
mass in the bulge volume of the RCGs measurement (±2.2 x ±1.4 x ±1.2kpc) with unprece¬ 
dented accuracy and robustness to be 1.84 ±0.07 x 10 10 M o . The stellar mass in this volume 
varies between 1.25 — 1.6 x 1O 1O M 0 , depending on the amount of dark matter in the bulge. 
We evaluate the mass-to-light and mass-to-clump ratios in the bulge and compare them to 
theoretical predictions from population synthesis models. We find a mass-to-light ratio in the 
K-band in the range 0.8 — 1.1. The models are consistent with a Kroupa or Chabrier IMF, but 
a Salpeter IMF is ruled out for stellar ages of lOGyr. To match predictions from the Zoccali 
IMF derived from the bulge stellar luminosity function requires 40% or ~ 0.7 x 1O 1O M 0 
dark matter in the bulge region. The BRAVA data together with the RCGs 3D density imply a 
low pattern speed for the Galactic B/P bulge of = 25 — 30kms _1 kpc -1 . This would place 
the Galaxy among the slow rotators ^ 1.5). Finally, we show that the Milky Way’s B/P 
bulge has an off-centred X structure, and that the stellar mass involved in the peanut shape 
accounts for at least 20% of the stellar mass of the bulge, significantly larger than previously 
thought. 

Key words: Galaxy: bulge - Galaxy: kinematics and dynamics - Galaxy: structure - Galaxy: 
center - methods: numerical 


1 INTRODUCTION 

Observations of external disc galaxies have shown that about half 
of all disc galaxies have strong bars ( [Eskridge et al.||2000| ). The 
Milky Way Galaxy (MW) has been considered for many years as 
one of them. The Galactic bar/bulge causes the non-circular mo¬ 
tions in the gas flow seen in Hi and CO (Peters, III|1975||Bmney] 
|et al.|1991[|Englmaier & Gcrhard| 1 999]|Fux| 1 999 \ and is the origin 

of the asymmetries seen in the near-infrared light distribution (Blitz 


& Spergel 1991; Weiland et al. 1994 

Binney, Gerhard & Spergel 

1997) and star counts (jNakada et a 

1991, Stanek et al. 1997; 


tic Bulge (GB) is regarded as the three-dimensional part of the bar 
seen nearly end-on ( |Shen et al.|2010[ |Martinez- Valpuesta & Ger-| 
|hard|201l| with a semi-major axis of about 2kpc (Gerhard 2002[ 
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[Wegg & Gerhard|2013| ), as is also indicated by the near-cylindrical 
rotation of the bulge stars ([Beaulieu et al.|2000||Kunder et al.|2012| 
|Ness et al.|20l3] ). 

In the last decade, stellar surveys of the GB such as 2 mASS 
( [Skrutskie et al.|2006) , VVV ( [Saito et al.|2012| ), OGLE ( |Sumi et aT] 
|2004| >, BRAVA ( [Rich et al.|2007| > and ARGOS ( [Freeman et al.|2013] > 
have released unprecedented data sets which allow us to study the 
GB star-by-star. The triaxial bulge of the MW is now believed to 
be a so called Box/Peanut bulge (B/P bulge) or X-shaped bulge as 
indicated by its bimodal distribution of Red Clump Giants (RCGs). 
This was first reported from the 2MASS catalogue by |McWilliam &| 
Zoccali p010|> an d from the OGLE-III survey by |Nataf et al. ( 2010| ). 


Ness et al.| ( |2012| > showed that this split red clump is seen for stars 


with metallicity [Fe/H] > —0.5. The peanut shape was mapped last 
year in three dimensions by | Wegg & Ge rhard ( 2013) using public 
data from the VVV survey. 

Star counts and infrared observations have revealed a long and 
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flat bar component, located mostly in the Galactic plane and ex- 
tendi ng up to / ^ 27° ([H ammersley et al. 2000, Benj amin et al.| 

|2005| |Cabrera-Lavers et al.|2007| >. Curiously, the angle of the bar 

relative to the line-of-sight to the Galactic Centre (GC) was in¬ 
ferred to be 0 = 45° in these studies, while the angle of the barred 
B/P bulge is accurately measured from the RCGs as (j) = 27° ± 2° 
jWegg & Gerhard|2013) . Study of N-body models suggested that 
the long bar and the B/P bulge data could be explained by a unique 
B/P bulge and bar structure formed by the buckling instability, if 
its two-dimensional outer bar component has developed leading 
ends through interaction with adjacent spiral arm heads ( |Martinez^| 
[Valpuesta & Gerhard|2011| ). 


The buckling instability is a well studied process in N-body 
simulations ( [Combes & Sanders|1981||Raha et al.|1991||Martinez-| 
Valpuesta & Shlosman 2004|). Cold stellar discs embedded in live 

dark haloes naturally tend to form a bar which experiences long 
term secular evolution because of angular momentum transfer from 
the disc to the halo ( [Debattista & Sellwood||2000| |Athanassoula| 
|2003) . During this secular evolution the bar can go through one or 
more buckling events ( |Raha et al.|1991 [[Martinez-Valpuesta, Shlos-| 
man & Heller 2006), making it vertically thick and creating the so 
called Box/Peanut shape, or X-shape in unsharp-masked images. 


In this context, the goal of this paper is to combine N-body 
simulations of evolved stellar discs with recent MW data to create 
dynamical models of the Galactic bulge and study its total mass, 
mass-to-light ratio, stellar mass and X-shape structure. To do so we 
use the 3D density of RCGs from [Wegg & Gerhard] ( |2013[ here¬ 
after also WG13) as well as kinematic data from the BRAVA survey 
( [Howard et al.|2008] |Kunder et al.|2012| ) to create particle models 
using the Made-to-Measure method (M2M). In M2M modelling, 
the weights of the particles in an N-body system are continuously 
updated until the observables from the model match a set of tar¬ 
get data constraints. The M2M method was introduced by [Syer &\ 
|Tremaine| \\996) and recast by |De Lorenzi et al.| ( |2007| > such that 
observational errors can be taken into account. We use the NMAGIC 
code o f|De Lorenzi et ak](|2007|) a dapted to the Milky Way prob¬ 
lem by Martinez-Valpuesta (2 0*12). NMAGI C has been used mostly 
in extragalactic studies (De Lorenzi et al.||2009[ |Das et al.~]|2011| 
Morg anti et al.|201~3) . M2M methods were previously used in the 


Milky Way context by |Bissantz, Debattista & G erhard ( 2004), Long 
|et al.| ( |20T3] ), also using BRAVA data, and |Hunt & Kawata| ( 201 4| ) . 


This paper is organized as follows. In Section [2] we present 
the set of N-body simulations of barred discs that we use as start¬ 
ing point for our M2M modelling. In Section [3] we recount the 
M2M method and the data sets we use to model the GB. The mod¬ 
elling is detailed in Section[4|where we present our best dynamical 
models and discuss the issue of the pattern speed of the bar. In Sec¬ 
tion [5] we show that we can recover the total mass in the bulge 
region with great accuracy, thereby relating the stellar mass to the 
dark matter mass in the bulge. In Section[6] we compute the stellar 
mass-to-light and mass-to-clump ratios of our models and compare 
them to theoretical predictions from population synthesis models. 
Section [ 7 ] quantifies the importance of the X-shape structure in the 
bulge which accounts for more than 20% of the stellar mass of the 
bulge. Finally, we discuss our results in Section [8] and summarize 
in Section [9] 


2 PARTICLE MODELS OF BARRED DISCS 

2.1 N-body models in different dark matter haloes 

Our M2M modelling of the GB relies on reasonable initial parti¬ 
cle models of barred discs. These initial models were created by 
evolving near-equilibrium stellar discs embedded in live dark mat¬ 
ter haloes. Near-equilibrium models are constructed using the pro¬ 
gram MAGALIE ( [Boily, Kroupa & Peiiarrubia Garrido||2001| ) and 
evolved with the tree-code GYRFALCON ( |Dehnen||2000| ),~~ail dis¬ 
tributed with the publicly available NEMO toolbox (Teuben 1995). 
During the evolution, the disc naturally forms a bar which rapidly 
buckles out of the Galactic plane and creates a B/P bulge ( [Combes| 
|& Sanders|1981||Raha et al.|1991| >. 

As we want to address the question of the total mass of the 
GB including the amount of dark matter in the bulge, we generated 
a set of five disc+halo N-body models using the same disc compo¬ 
nent and varying the halo properties. The disc is exponential with 
scale length of 1 internal units (iu), scale height of 0.14iu, unit to¬ 
tal disc mass, and Q parameters of 1.4 at R = 3.07iu. Halos have 
a Hernquist density profile with flattening of 0.8 and a sharp cutoff 
at 20 iu. All models contain two million particles, one million each 
for disc and halo. With these settings fixed only two free parame¬ 
ters remain: the halo mass inside the cutoff, M^, and the halo break 
radius of the Hernquist profile, a^. They were fixed considering the 
total rotation curve of a model. Following the language of |Sackett| 
\\991) , we call the degree ofmaximality of a disc the proportion of 
the disc contribution to the total velocity curve, at the radius where 
the disc velocity curve is maximal. Assuming a flat rotation curve at 
large radii, we build a one parameter family of models parametrized 
by the degree of maximality of the disc. We make five models of 
this family, called M80, M82.5, M85, M87.5 and M90 which have 
different degree of maximality ranging from 80% to 90%. As we 
show later this allows us to build models of barred galaxy which 
span all kind of rotators, from slow to fast. The halo parameters 
used to construct these models are summarized in Table [T] 

For each model, we select a snapshot of a late evolutionary 
stage, where the bar is fully grown. The circular velocity curves ob¬ 
tained from the azimuthally averaged potential of this set of models 
are plotted in Fig. [I] before any evolution (top) and, for the selected 
snapshot after bar and B/P bulge formation (bottom). The different 
colours refer to different models as stated in the legend. In the up¬ 
per plot, the solid black line is the circular velocity curve of the 
initial disc component (common for all models) and the vertical 
dashed line indicates the radius at which the degree of maximality 
is determined. As expected the bar formation moves material in¬ 
ward which increases the circular velocity in the inner region of the 
disc. 


2.2 Geometry and scaling 

According to the latest studies ( [Chatzopoulos et al.|2015[|Reid et al.| 
[2014| ), we assume a distance to the Galactic Centre (GC) of Rq = 

8.3 kpc. The bar is placed at an angle of a = 27° from the Sun-GC 
line (WG13). All our models are scaled independently based on 
the length of their long bar. Several studies based on different data 
sets indicate that the long bar of the MW ends at about / = 27°. 
Hammersley et al. ( 2000) studied star counts of K giants in several 
in-plane fields and places the end of the long bar at / ~ 27° while 
|Cabrera-Lavers et ak| ( |2007[|2008) found respectively values of 27° 
and 28°, based on RCGs counts from the 2 mASS catalogue and 
from UKIDSS data. |Benjamin et ak| ( |2005| ) also found indication of a 
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Figure 1. Top : Circular velocity curves of our five disc-halo models before 
evolution, in model units. The black line shows the circular velocity curve 
of the disc only, kept the same for all models. The dashed vertical line in¬ 
dicates the radius where the maximality of the disc is computed. Bottom : 
Azimuthally averaged circular velocity curves after the bar and B/P bulge 
formation, also in model units. 


long bar for / < 30° using star counts from the GLIMPSE catalogue. 
All these different estimates agree well with each other and give 
therefore a consistent length scale. With Rq = 8.3kpc, a bar angle 
of 27° and the end of the bar at l = 27°, the half length of the bar is 
/?bar — 4.66kpc. The different parameters quoted above constitute 
our fiducial set of parameters. We checked that our main results are 
not significantly affected by assuming a bar angle of 32° instead of 
27° or a shorter bar half-length of /?bar = 3.8kpc (see Section [5^2| , 
or by using Rq = 8 kpc instead of 8.3 kpc. 

We used ellipse fitting of the face-on projection of our models 
to compute their bar lengths, and scaled them to 4.66kpc. The ve¬ 
locity scaling is kept free and will be determined dynamically from 
the data during the modelling process, as explained in Section [34] 

In order to avoid referring to scale dependent quantities we 
will refer to the rotational velocity of the bar using the dimen¬ 
sionless number = R C r/Rb ar> the ratio between the corotation 
radius and the half-length of the bar. Bars with & > 1.4 are called 
slow rotators, while those with ^ ^ 1.4 are fast rotators ( [Debattista] 
|& Se llwood 20 00]). Beca use the bar cannot extend beyond corota¬ 
tion ( Contopoulos|1980| , & has to be larger than one. Our models 
span & values quite uniformly from 1.8 to nearly 1.1, which cor¬ 
responds to the full range of reasonable values for barred galaxies 
JEhnegreen 1996). Statistics for external galaxies from[Rautiainen, 
Salo & Laurikainen ( 2008) show that nearly all galaxies of Hubble 
types SO, SB a or SBb are consistent with being fast rotators while 
for SBc galaxies (like the MW), bars can be either fast or slow. The 
sampling of & values we obtain is a consequence of the different 
haloes we used. The halo absorbs angular momentum during the 
bar formation and therefore the more maximal the initial disc, the 
less halo material there is to absorb angular momentum and so the 
faster the bar. One should keep in mind that for our set of models, 


Table 1. Main parameters of our five disc-bar-halo models. The two first 
rows show the initial halo parameters: the halo mass inside cutoff in 
units of the disc mass M d and the Hernquist break radius a^. The next rows 
give parameters after bar and bulge formation and scaling to physical units 
(see Section |L2) : the corotation radius R cr of the bar, the ratio of the corota¬ 
tion radius over bar half-length and the dark matter fraction in the bulge 
MdmMoi. 



M80 

M82.5 

M85 

M87.5 

M90 

M h /M d 

9.42 

8.51 

7.68 

6.94 

6.27 

a h [m\ 

17.62 

19.44 

22.05 

26.16 

33.68 

Rcr [kpc] 

8.4 

7.6 

6.7 

5.6 

5.0 

& Rcr/ Rbar 

1.80 

1.64 

1.44 

1.21 

1.08 


44.2% 

40.8% 

34.5% 

28.6% 

25.2% 


the halo mass in the inner parts and the & value are not indepen¬ 
dent parameters. The different model and bar parameters are given 
in Table [T] 

Throughout this paper the (x,y,z) frame refers to the Galac- 
tocentric inertial frame where z is the vertical axis pointing to the 
Galactic North and y the Sun-GC axis. The bar rotates at pattern 
speed £2 p in this inertial frame and we shall refer to the rotating 
frame of the bar as (x\y\z f ), with x',/ and z! respectively the ma¬ 
jor, intermediate and vertical axis of the bar. Fig.[2]shows the face- 
on and side-on projections of our models with these geometry and 
scaling. 


3 MADE-TO-MEASURE WITH MW OBSERVABLES 


As an alternative to distribution function-based methods (Dejonghe 


1984] |Binney||20 


Illingworth 1990 


0]), moment-based methods fBinney, Davies & 


_ fCappellari et al.||2009| or orbit-based meth¬ 

ods (|Schwarzschild|1979||Thomas et al.|2009) , |Syer & Tremaine| 
(1996]) proposed a particle-based algorithm to study stellar dynami¬ 
cal equilibria, known as the Made-to-Measure (M2M) method. This 
algorithm consists of adapting the particles weights of an initial par¬ 
ticle model of the system of interest such as to reproduce a given 
set of observables. |De Lorenzi et al. ( |2007| hereafter DL07) im¬ 
proved the original method by Syer & Tremaine (1996) to take into 
account observational errors, and implemented it as the NMAGIC 
code. NMAGIC has been used in numerous studies, mostly in the 
context of elliptical galaxies (e.g. |De Lorenzi et al.|2009[|bas et al.| 
2011; Morganti et al. 2013) and has been adapted to barred galaxies 
by Martinez-Valpuesta (2012). In the context of the Milky Way, the 
M2M method has previously been used by Bissantz, Debattista &| 
|Gerhard| ( [2004] ), |Long et al.| ( [20T3] ) and |Hunt & Kawata| ( |2014| ).~ 


3.1 Theory of the M2M method 

Let us consider a system characterized by its DF f(z) defined on 
the phase space of the system. Any observable yj of this system 
will be written as 

yj = J Kj(z)f(z)d 6 z, (i) 

where Kj is the kernel of the observable and z the phase space vec¬ 
tor. If we represent the system by a set of N particles, with particle 
weights Wi(t), the observable will be written as 

yj(t) = 'LK j (z i (t))w i (t), ( 2 ) 

i =1 
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Figure 2. Face-on and side-on views of our five initial models of barred 
discs with B/P bulges. The Sun is located 8.3kpc from the Galactic center 
and the bar is at an angle of 27° from the Sun-GC line. The dotted lines on 
the face-on view indicate sight lines spaced every 10° in Galactic longitude. 
The bold rectangle indicates the boundary of the box where the density is 
given by Wegg & Gerhard|(2013) , see Section [±2| 


where Zi(t) is the phase-space coordinate of particle i at time t. The 
Wi are proportional to the physical weights of the particles but one 
can also see them as density-elements of the phase-space. 

Let us now consider different data sets, indexed by the sub¬ 
script k, that one wants to fit. The difference between the model 
and the observational target is quantified using the residuals A k {t) 
defined as 


A k j(t) = 


y^hll 

<*{Yj ) 


(3) 


where Y k denotes an observable of data set k, c(Yj) its associated 
error, andy^(f) the corresponding model observable. 


The M2M method will adapt the weights of the particles in 
order to maximize the profit function F defined by 


F = ^S- l -xl v 

where 


i 

and 


^t = l4(A k j) 2 . 

kj 


(4) 

(5) 

(6) 


S corresponds to an entropy term used to regularize the particle 
weights in order to ensure that they do not deviate too much from 
a set of predefined priors w*. The are numerical weights of the 
different data sets, as formally introduced by |Long & Mao ( 2010) ). 
The determination of these factors is discussed in Section U.6.31 
The heart of the M2M method is the following weight evolu¬ 
tion equation: 


dwi 

dt 


= ewi(t ) 


dF_' 

dwi 


= emit) 


** - E E }} a j (o 

dm t j CT (Yf) 


(7) 

( 8 ) 


where the bracket term is the so-called force-of-change. Note that 
the passage from Eq.[7]to Eq.[8]is made under the assumption that 
the kernels K k - do not depend on the weights of the particles. We 
will make this assumption as it allowed Syer & Tremaine (1996) 
and DL07 to prove the convergence of the particle weights for small 
linear deviations from the solution. No additional term is used in 
Eq.0 in particular we do not re-normalize the weights of the par¬ 
ticles. This is done on purpose to allow fitting different masses of 
the bulge as shown in Section [4] 

In order to reduce the shot noise of the particle model we fol¬ 
low [Syemfe^emaine (1996) and DL07 and artificially increase the 
effective number of particles using temporal smoothing, replacing 
yj(t) in Eq.[3]by 

yji*) = Jyj{t-T)e~ aT dT . (9) 


3.2 Density observables 

3.2.1 3D density of Red Clump Giants 

Using the Red Clump Giants (RCGs) from the VVV survey, WG13 
measured the three-dimensional density distribution of the inner 
part of the Galactic bar/bulge. They evaluate line-of-sight den¬ 
sity distributions of the RCGs by deconvolution of extinction and 
completeness corrected K s band magnitude distributions for dif¬ 
ferent VVV fields. To do so, they fit a background to each K s 
band magnitude distribution and identify the RCGs as the ex¬ 
cess over this background. Assuming an 8-fold mirror symme¬ 
try, they constructed a 3D density map of the RCGs in the inner 
±2.2 x ±1.4 x ±1.2kpc along the (x',y',z') axes of the galactic 
bar/bulge, and checked that the data show only small deviations 
from 8-fold symmetry (their Fig. 15). Their work relies on several 
assumptions that they carefully investigated for systematic varia¬ 
tions of their fiducial parameters. They finally provide us with one 
fiducial and five variant density maps of RCGs in the bulge. We use 
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the range of these variant maps as a systematic error on the density 
measurements. The typical magnitude of these errors is about 10%. 

We assume here that RCGs are good tracers of the stellar mass 
in the bulge, i.e. that the ratio of number of RCGs per unit of stellar 
mass stays constant in the bulge. This assumption is supported by 
the theoretical work of |Salaris & Girardi] ( |2002| ) where they showed 
that for an old star population of about lOGyr one would expect 
the number of RCGs to vary by less that 10% for metallicity in 
the range —1.5 < [M/H] ^ 0.2. As the bulge appears uniformly 
old with no significant metallicity component extending beyond the 
range —1.5 ^ [M/H] ^0.2 (Zoccali et ak]2003]>, we can consider 
RCGs as good tracers of the stellar mass and therefore use the map 
of WG13 as a constraint on the shape of the stellar mass density in 
our models. Note that the map gives only the shape of the density 
and not its absolute value because the number of RCGs per unit of 
stellar mass is a priori unknown. 

Unfortunately high extinction and crowding prevented WG13 
to reliably measure the RCGs density within 150pc of the Galac¬ 
tic plane, and also cause some uncertainty immediately above this 
150pc strip. Using directly the original map with this missing strip 
would be inappropriate since we want to model different bulge 
masses by changing the scaling of the density constraint. An in¬ 
complete map would let the midplane free and lead to unrealistic 
models, with for example a light in-plane component and a massive 
out-of-plane component. Hence we model the density in the miss¬ 
ing ±150pc strip by extrapolation of each vertical density profile 
of the original 3D map. 

We found that the vertical density profiles in both the data 
of WG13 and the initial particle models are well represented by a 
sech 2 function. Our fiducial extrapolation is thus based on the best 
sech 2 fit of each vertical profile. We account for the uncertainty due 
to the choice of the extrapolation law by considering a different one 
in Section |5.2[ showing that the total bulge mass is insensitive to 
this choice. 

The fiducial extrapolated density map of RCGs in the ±2.2 x 
±1.4 x ±1.2kpc box is plotted in projection in Fig. [3] The side-on 
view shows a very strong peanut-shape. 


3.2.2 Implementation in NMAGIC 


The 3D density map was evaluated on a regular Cartesian grid of 
(30,28,32) cells along the (x / ,y / ,z / ) axis, covering a box of the in¬ 
ner ±2.2 x ±1.4 x ±1.2kpc, as in WG13. We will refer to this re¬ 
gion as the “bulge-in-box” or the abbreviation b-b. From WG13 we 
know the total number of RCGs in the b-b, but the corresponding 
stellar mass must still be determined from dynamical modelling. 
We parametrize the total stellar mass of the b-b using a dimension¬ 
less factor, defined as the ratio of the target stellar mass of the 
b-b divided by its value in the initial model. The target density ob¬ 
servables Yj, corresponding to the target stellar mass in cells j, in 
internal units, are then given by 




d - ^ I £ Wl (t = 0) 
ie b-b 


WRCGsOQA 3 * 

/b-b n RCGs d 3 X 


( 10 ) 


where the sum is over all initial particle weights in the b-b, wrcgs in 
the number density of RCGs, and A 3 x the volume of the cells. For a 
given scaling from model internal units to physical units, a change 
in the value of & is equivalent to a change in the target stellar mass 
of the b-b. For a given NMAGIC takes care of increasing or de¬ 
creasing the weights of the particles to reach the target observables 
Yj, thereby modelling different-mass bulges. 


£«c[pc] 2 

0.1 0.3 0.9 2.9 8.8 



Figure 3. End-on (top left), face-on (top right) and side-on (bottom right) 
projection of the extrapolated 3D RCGs density map originally from |Weg~g| 
|& Gerhard](2013) . The dashed lines in the end-on and side-on view show 
the ±150 pc region where the density map was extrapolated. 


The model density observables in the cells of the b-b are then 
computed from Eq. fusing the following kernel 


Kj(Zi) 


1 if i e cell j, 
0 otherwise. 


( 11 ) 


3.3 Kinematic observables 

3.3.1 BRAVA data 

We kinematically constrain our models using data from the Bulge 
RAdial Velocity Assay (BRAVA) { [Rich et aL]|2007| |Howard et al.| 
2008 , |Kunder et al.|2012| >. The BRAVA survey is a large spectro¬ 
scopic survey of M giant stars selected from the 2MASS catalogue. 
According to Howard et al. ( 2008), the light of these M giants traces 
the 2qm light of the GB and therefore M giants are a good proxy to 
study the overall kinematics of the GB. The BRAVA survey provides 
us with galactocentric rest-frame mean radial velocity and velocity 
dispersion in more than 80 fields through the bulge, mostly between 
/ = —10° to / = 10°. On starting this project the ARGOS data ( |Ness| 
|et al.|2013] > were not yet available so we restricted ourselves to the 
BRAVA data. The ARGOS data have a more complicated selection 
function and will be included in a later paper. 

Before constraining a model it is important to study in detail 
how the BRAVA stars were selected in order to reproduce any se¬ 
lection bias. The BRAVA stars were selected from their location on 
a K versus J — K colour-magnitude diagram aiming to select only 
bulge members with no metallicity bias. |Howard et al.[ ( |2008|) de¬ 
ployed a lot of effort to ensure no metallicity bias but the bulge 
membership criteria based on magnitude cuts adjusted by eye are 
more questionable. We used the GALAXIA model ( [Sharma et al.| 
[201 1| ) to evaluate the possible foreground contamination using the 
same selection criteria. We found that contamination was negligi¬ 
ble toward the centre but appears to rise with increasing Galactic 
longitude, reaching 20% of the sample at / ~ 15°. We found no sig¬ 
nificant variation with latitude in the latitude range of the BRAVA 
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fields. We decided not to simulate this contamination, given the fact 
that the foreground discs in our models are not made to match the 
disc of the Milky Way. Instead, we simply exclude from our analy¬ 
sis all the fields outside the inner ±10° in longitude, where contam¬ 
ination is probably significant. The magnitude cuts also introduce 
a slight bias towards the near side of the bulge. Faint M giants are 
more numerous than bright ones so have a larger probability to be 
part of the sample on the near side than on the far side of the GB. 
We model this effect below. 


3.3.2 Implementation on NMAGIC 

As stated in SectionJTT] the linear convergence of the M2M method 
is assured if our observables are of the form given by Eq.[2] where 
the kernels do not depend on the weights of the particles. Therefore 
as our observables, instead of using the mean radial velocity and 
velocity dispersion we use the first and second mass weighted radial 
velocity moments, indexed as v r i and v r 2 . The observables are then 
of the form of Eq.[2]with the following kernels: 


a check of our modelling, comparing the data to our model pre¬ 
dictions, without fitting them. Rattenbu ry et al.| < |2007a| > selected 
bulge RCGs from their colours and magnitudes, and excluded all 
stars with total proper motion larger than lOmasyr -1 . The bulge 
membership of model particles is evaluated using the same crite¬ 
rion as for the BRAVA observables, i.e. |y| ^ 3.5kpc. Similarly to 
the BRAVA observables, our proper motion observables are the first 
and second velocity moments in the / and b direction, indexed by 
v/i 2 and v^i 2 - The kernels of the first and second velocity mo¬ 
ments in the / and b directions are similar to those of Eqs.[T2|and 
13 replacing v-" by the velocity of particle i in the / and b direction, 
expressed in mas yr -1 . 

Our selection function when evaluating the model proper mo¬ 
tion observables is then given by: 


d] l ^\ Zi ) = < 


1 if i G field j, |y*| < 3.5kpc 
and ^v* + vg ^ lOmasyr -1 
0 otherwise 


(15) 


Ky i (z i ) = 8] r '(z i y i d2) 

Kp( Zi ) = 8p(zi)(^) 2 (13) 

where v •" is the radial velocity of the particle i and <5j H ,Vr2 are the 
field selection functions. In order to remove foreground contam¬ 
ination we consider only particles whose y coordinate (along the 
GC-Sun axis) is in absolute value lower than 3.5 kpc. This corre¬ 
sponds to the selection criterion used in the ARGOS survey ( |Ness| 
|et al.|[2013| and its use here is supported by the fact that ARGOS 
data and BRAVA data agree with each other. We checked that the 
exact form of this selection function does not change significantly 
the kinematic observables in the inner 10°. 

In order to map the bias toward the near side we use the ap¬ 
proximate luminosity function of giant stars in the bulge d>(M^) oc 
IqO.28 m k f rom WG13 where Mk is the absolute magnitude in the 
K band. In each field, stars of the BRAVA sample are uniformly se¬ 
lected between two magnitude cuts from their apparent magnitude 
mx = Mk + 51og(r/10 pc) + Ak where r is the distance to the Sun 
and Ak is the extinction. If we consider the extinction as a fore¬ 
ground extinction, this uniform selection in is equivalent to a 
non-uniform selection in r along the line of sight with weighting 
1Q-O.28x5log(r) _ r -l.4 This r -1 ' 4 weighting lowers the natural r 2 
weighting due to the cone opening of the line-of-sight. We therefore 
adopt the following selection function: 


\zi) 


r- 1-4 if i G field j and |y;| < 3.5kpc, 
0 otherwise. 


(14) 


In order to compare the model observables to the target data, 
we have to weigh the target data by the expected mass in each field. 
As the stellar mass in each field is unknown, we use the model 
mass, and update this weighting several times during the fit (see 
Section |T6] >. 


3.3.3 Proper motions 

In addition to the BRAVA data, we use proper motion data from 
Rattenbury et al. (2007a). These authors computed proper motion 
dispersions, a/in the / and b directions for a large number of 
bulge RCGs in 45 bulge fields from proper motion measurements 
for stars in the OGLE-II survey. We use these proper motions as 


3.4 Dynamical velocity scaling 


Our models are evolved in a system of internal units where the 
length unit, velocity unit and gravitational constant are set to unity. 
When comparing model to data we scale the models to physical 
units using the length of the bar (see Section [22] ), the gravitational 
constant and a velocity scaling. This velocity scaling is first fixed 
to some referen ce value using t he circular velocity at the radius of 
the Sun from Bovy et al. (2012) (218 kms -1 ). Keeping the velocity 


scaling fixed to this value would be inappropriate, however: our 
models are constrained in the inner part by the BRAVA data and the 
3D density but no effort has been made to make the model rotation 
curves match the MW rotation curve at larger radii. Therefore the 
velocity scaling given by the circular velocity at the radius of the 
Sun is not relevant in the bulge. Hence we determine the velocity 
scaling directly from the BRAVA data using a variant of the method 
presented by |De Lorenzi et al.] ( |2008| ). In Eqs.[l2]and[T3]we replace 
the radial velocity of a particle v[ by yv-", where v \ is now the radial 
velocity of particle i expressed in physical units using the reference 
velocity scaling, and y is a numerical factor initially set to one. 
Through Eqs.gjandjijthe total % 2 (Eq. [6j and the profit function F 
(Eq. 0 now depend on y. To find the maximum of F with respect 
to y, we use the following evolution equation for y, similar to the 
force-of-change equation: 


dy d% 2 


(16) 


where p sets the magnitude of the force-of-change applying on y. 
The velocity scaling plays a role for the kinematic observables only, 
so the derivative of X 2 is given by 


^=2£A, 
dy u 


-V 1 


dAY 

~W +K 


dA ]' 2 

A ry 


(17) 


where the derivatives of the A I 1,2 are given by the following equa¬ 
tions: 

. 


,, v rl 


dy yo(Yj rl ) 


d _Al 

dy 


Vrn 

2 >f 

V r 2\ 


r<r(r, 


(18) 

(19) 
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The value of 17 should be large enough to ensure the conver¬ 
gence of / during the fit, but the time-scale for its evolution must be 
longer than the temporal smoothing time-scale. We fixed 17 to 0.4 
from estimating the magnitude of the right hand side of Eq. |T6] In 
all our fits, 7 converges to some final value which we use to convert 
internal units to physical units. All velocities and proper motions 
are scaled by 7 with respect to their value with the reference scal¬ 
ing and all masses are scaled by y 2 . 


3.5 Potential and model integration 

Stellar and dark matter particles are all integrated in their com¬ 
bined potential that comes initially from our models of barred discs 
evolved in dark matter haloes (see Section[2]). Only stellar particles 
are used as M2M particles, whose masses are both used as gravi¬ 
tational masses and as M2M weights. During the NMAGIC fits, the 
gravitational potential changes as the particle weights in the bulge 
adapt to match the target density (including the & factor) and the 
kinematic constraints. To take this into account, the potential of the 
stellar and halo particles is recomputed from time to time during 
the fit. This ensures that the new weights of the disc particles can 
influence the kinematics of the model and that the final converged 
model evolves in its own gravity. In between potential updates, the 
particles are integrated in a frozen rotating potential. The potential 
always rotates at the constant pattern speed £ 2 p of the initial model, 
in internal units, and the centre of mass as well as the rotation axis 
of the potential are kept fixed during the complete NMAGIC fit. Each 
potential update corresponds to an update of the orbit library. 

In all our fits the potential is recomputed 10 times during the 
weight evolution, which we found is sufficient to ensure the self¬ 
gravity of the converged model as well as smooth updates of the 
orbit library. The potential is computed using the 3D polar grid 
code from Sellwood & Yalluri (1997). The particles are integrated 
with a drift-kick-drift adaptive leap-frog algorithm. The integration 
scheme is such that over a typical fit integration time in a fixed 
potential, the Jacobi energy is conserved to a level of 10 -3 or better. 


3.6 NMAGIC parametrization 

3.6.1 Fitting procedure 

A typical NMAGIC fit consists of the three following phases. First 
we evolve the particles T smoot h iterations during which we com¬ 
pute the model observables in order to initialize the temporally 
smoothed observables. Then we integrate the model for Tm 2 M it¬ 
erations, changing the weights of the particles according to Eq. [ 8 ] 
while updating the potential regularly. All observables are matched 
to the data at the same time. After a last potential update we finally 
relax the model for T re i ax iterations without changing the particle 
weights. This last phase is important to avoid over-fitting the data 
and obtain realistic models. The usually slight % 2 increase during 
this relaxation reveals how much the model was forced to fit the 
data by Eq. [ 8 ] At the end of the run, we also checked the con¬ 
vergence of the particle weights, using the convergence criterion 
detailed in |Long & Mao| ( [2010| ). A particle weight is considered 
to have converged if its maximum relative deviation from its mean 
value over some period of time is smaller than some threshold. As¬ 
suming a period of time corresponding to four circular orbits at 
2kpc and a threshold of 10%, about 97% of the particles weights 
converge in a typical NMAGIC fit. 


Table 2. Typical set of NMAGIC parameters. T smoot h, Tm 2M and T re i ax are 
the number of iterations for the 3 phases of a NMAGIC fit. T$ and T mass are 
the number of iterations between two potential updates and updates of the 
mass weighting for the kinematic observables. 1 /(a*dt) is the time-scale 
of the temporal smoothing in number of iterations, e/wo is the magnitude 
of the force of change and wo is the initial weight of the stellar particles. 
The A parameters are described in Section [T.6.3| 


Tsmooth [it] 

Tm 2M [it] 

Treiax [it] 

T(j> [it] T mass [it] 

10 4 

10 5 

2x 10 4 

O 

O 

1 /(a*dt) [it] 

e/w 0 [ iu 1 ] 

'td 

^Vrl,V r2 

2.5 x 10 3 

0.04 

1 

25 


3.6.2 Parameter values and time-scales 

The parameter values we use are shown in Table [2] The iteration 
step of the M2M procedure dt is fixed to one thousandth of the 
time needed to complete a circular orbit of radius 2kpc. All mod¬ 
els are integrated for a constant number of dt so the number of 
bar rotations during the weight adaptation phase varies from model 
to model, ranging between 25 for model M80 and 35 for model 
M90. This corresponds to a physical integration time between 6 
and 7 Gyr, once scaled to physical units using the velocity scaling 
determined dynamically by NMAGIC. l/(a * dt) is the time-scale 
of the temporal smoothing expressed in terms of number of M2M 
iterations, and e/wq is the magnitude of the force-of-change in Eq. 
[ 8 ] where wo is the initial weight of a stellar particle. 

3.6.3 Xfr parameters and regularization 

The different set of observables, here RCGs 3D density and BRAVA 
data are weighted by the A^ parameters in the profit function. The 
force-of-change (Eq. [ 8 j already takes into account observational 
errors so in theory these A^ should all be set to 1 , in order to re¬ 
ally minimize the total % 2 . In practice, experiments showed that 
with all Afc = 1 the model ignores completely the kinematic con¬ 
straints. This is mostly caused by the very large number of den¬ 
sity constraints (26880) with respect to kinematic constraints (164). 
As we do want to fit the BRAVA data, we increase the weight¬ 
ing of BRAVA constraints with respect to the RCGs density using 
Ay, / Ad = A Vr2 / Ad = 25, even if then we no longer strictly minimize 
the total x 2 - These A^ were determined using the distribution of the 
force-of-change contribution of each set of observables, such that 
the mean force-of-change due to the density observables should be 
equal to the mean force-of-change due to the BRAVA observables. 
This causes the BRAVA data to be reasonably fitted without being 
over-fitted. A strong over-fitting would lead to a significant increase 
of x 2 during the relaxation step at the end of the fit. 

We found that our models were smooth enough without using 
any entropy smoothing. This is a consequence of the very dense 
constraint from the 3D density data. Hence we chose to set p, = 0 
in Eq.[ 8 ]for all our models. 

4 DYNAMICAL MODELS OF THE MW 

In this section we fit our initial models to the BRAVA data and the 
3D RCGs density using NMAGIC. The models differ by their dark 
matter fraction in the inner part and by their dimensionless corota¬ 
tion radius & (see Section [O] ). In Section |4J~1 we determine the 
best bulge stellar mass & for each model, using the full NMAGIC 
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modelling procedure described above, and in Section pk2| we com¬ 
pare all models with their respective best &. In Section |4.3[ we 
constrain the pattern speed of the MW bar from this modelling, and 
in Section |4~4| we compare the models with available proper motion 
data. 


4.1 Determination of the stellar mass in the bulge 


In Section |372| we parametrized the target density observables using 
a free numerical factor &. The value of & directly sets the stellar 
mass of the bulge in internal units for each model. We find that & 
has a strong influence on the velocity dispersions but only a very 
slight influence on the mean velocities. Empirically, the shape of 
the mean velocity profile in our models is fixed by the density dis¬ 
tribution. The velocity amplitude in internal units is then fixed by 
the pattern speed (Eq. 8 of |Debattista, Gerhard & Sevenster|2002| . 
The larger kinetic energy in models with larger bulge mass & can 
therefore only be put into the velocity dispersion. This is shown 
in Fig. [4] where the kinematic observables of model M85 are plot¬ 
ted along with the BRAVA data. In this figure the different coloured 
lines show the kinematic profiles obtained after fit of the density 
normalized with different values of &. The upper (lower) plots 
show mean radial velocity (velocity dispersion) profiles for three 
different latitudes as well as along the minor axis. For the compar¬ 
ison, the model observables in Fig. [4] have all been plotted for the 
same scaling constant so as to match only the mean velocity data. 
As & has nearly no influence on the velocity, this better highlights 
the effect of & on the dispersion. 

Fig. [4] illustrates how more massive bulges lead to higher ra¬ 
dial velocity dispersions. By finding the value of & which gives 
the best agreement with the data, we can recover the stellar mass of 
the b-b for each model. Fig. [ 5 ] shows the y} per data point of the 
velocity and velocity dispersion for all models plotted in Fig.[4] ver¬ 
sus &. As expected, a clear minimum is present in the dispersion 
plot and provides us with a best value of & for each model. Best 
values of & are given for all models in the first column of Table 
[3] With the optimal velocity scaling determined in these NMAGIC 
fits, these & values can be converted into physical values of stellar 
mass in the b-b. We find values of 1.25 x 1O 1O M 0 for model M80 


IVlQ IQ 

ion [ 5 ] 


to 1.6 x lO iU M 0 for model M90; see also Section 

& was sampled on a regular grid with spacing 0.1 and is there¬ 
fore determined with an accuracy of 0.05. Due to the rescaling dur¬ 
ing the NMAGIC fits, an uncertainty of 0.05 in & typically corre¬ 
sponds to a change of less than 0.025 x 1O 1O M 0 in the final stellar 
mass and less than 0.01 x 1O 1O M 0 in total mass (stellar and dark 
matter in the b-b). 


4.2 Best dynamical models of the Milky Way bulge 

Now we compare our five best dynamical models M80-M90 with 
different dark matter haloes, obtained after NMAGIC fit to the data 
for their respective best value of &. For all models the density 
and its peanut shape is well fitted as shown in Fig. [6] This figure 
compares the contours of the three principle axis projections of the 
density in the bulge after fitting with the target RCGs density. Each 
colour line represents one of our best mass models and the black 
line is the projection of the target density. 

The radial velocity and dispersion profiles of the same best 
& models are plotted in Fig. [7] Except for small differences in 
the shapes of the dispersion curves, they all look very similar and 
provide a good fit to the data. This similarity is explained by the fact 



0.6 0.8 


M80 


1.0 


1.2 

& 


1.4 


1.6 


1.8 2.0 


IZZI M87.5 


H M82.5 IZZI M85 
■ M90 


Figure 5. % 2 per data point of the velocity observables (upper plot) and 
velocity dispersion observables (lower plot ) as a function of the bulge mass 
factor &. Each colour represents a different initial model from Section[2]as 
shown in the legend. In all cases, a best value of & is clearly visible from 
the velocity dispersion x 2 - -curves. 
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Figure 6. Contour plot of the projected 3D density in the bulge of our 
best dynamical mass models compared to the measured RCGs density from 
|Wegg & Gerhardj {2013) . The projections are the same as in Fig.[3]and the 
contours are spaced by a third of a magnitude. In all cases the density is 
well fitted. 


that the shape of the velocity profile is mostly fixed by the shape of 
the density and its magnitude is adapted to the data by the floating 
velocity scaling. The magnitude of the velocity dispersion profiles 
can then be adapted independently by the & parameter. 

More quantitatively, the % rl per data point is shown in Table [ 3 ] 
separately for the density, the total kinematics, the velocity only and 
the velocity dispersion only. Values in this table are not weighted by 
the corresponding parameters. The total /^-weighted chi square 
Xlt/ntot = E ^kXk /H^k n k actually minimized by NMAGIC is given 
in the last column of Tabled 

All our models provide good fits of the RCGs density with a 
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Figure 4. Comparison of the brava data to the fitted model M85, for different values of the stellar mass parameter The upper (lower) plots show the 
BRAVA mean radial velocity (velocity dispersion) profiles for b = —4°, b = —6°, b = — 8° and along the minor axis / = 0°. The different colours indicate 
different values of & as stated in the legend. The influence of & on the kinematics is clearly visible in the dispersion plots. The small inserts on the top right 
of each plot show the BRAVA fields in Galactic coordinates and highlight the fields shown in the corresponding plot. 



Figure 7. Velocity and velocity dispersion profiles for our five best dynamical models of the Milky Way with different dark matter haloes. The plotting 
conventions are the same as in Fig.[4]except that here the colours indicate different models as stated in the legend. 
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Table 3. Best values of & and % 2 per data point of the density observables 
(Xd/ n d)* t° ta l kinematic observables (X^/ n k), velocity only (Xv/ n v) and 
dispersion only (% 2 /%) for our five best dynamical mass models. The last 
column shows the reduced weighted % 2 actually minimized by NMAGIC. 


Model 

& 

Xi/ni 

*k/«k 

Zv/«V 

ill no 

Xlt/ntot 

M80 

1.0 

0.42 

1.45 

1.67 

1.24 

0.56 

M82.5 

1.1 

0.32 

1.46 

1.68 

1.24 

0.47 

M85 

1.1 

0.24 

1.46 

1.65 

1.27 

0.40 

M87.5 

1.3 

0.23 

1.52 

1.73 

1.31 

0.40 

M90 

1.5 

0.30 

1.52 

1.75 

1.29 

0.46 


final Xd/ n d from 0.23 to 0.42. As the errors in the density are sys¬ 
tematic, one should not over-interpret these values. The kine¬ 
matics are also well fitted for all models with the model disper¬ 
sions marginally steeper than the data for large |/| and \b\, and 
Xk/ n k ranging from 1.45 for model M80 to 1.52 for model M90. 
This kinematic x 2 is significantly better than that obtained by |Long| 
[et al.||20l3] ). Given the scatter visible in the kinematic data, these 
models are good candidates to represent the MW bulge even if the 
kinematic x£ is about 1.5. 

Because the different models all give a reasonable fit, we do 
not rule out some of them based on simple x 2 considerations. 
One has to be careful when drawing conclusions comparing x 2 
values on typical M2M problems. Indeed, as shown by [Morganti| 


et al. 


2013b the common practice to use Ax 2 = X 2 ~ Xmin an ^ X 2 


statistics to evaluate confidence levels on x£ 


min 

is not appropriate 

for their M2M results. The Ax 2 analysis requires a positive num¬ 
ber of degrees of freedom, while it is usually negative in M2M 
problems because the particles are vastly more numerous than the 
data constraints. Hence, the X 2 ! n given in Tableware only an 
indication of the distance between model and data. These values 
are also strongly influenced by a few data points which appear to 
be possible outliers. This is the case for example in the fields at 
(l,b) = (—7°, —6°); (—5°, —4°); (—4°, —8°). Removing these pos¬ 
sible outliers reduces the absolute value of the kinematic X 2 / n k by 
about 0.2. It does not however change our best value of &. 

The face-on and side-on projections of our final best mass 
models are plotted Fig. [8] Even though not enforced by NMAGIC, 
the long bar component is still there in the fitted models. Its pres¬ 
ence indicates that the gravitational potential updates performed 
during the fit were smooth enough to keep long bar particles on 
bar orbits. 


4.3 The pattern speed 

The pattern speeds of our models are converted to physical units 
using the velocity scaling determined by NMAGIC in order to fit the 
BRAVA data. For models M80, M82.5, M85, M87.5 and M90 we get 
the values £2 p [kms -1 kpc -1 ] = 24.7, 25.7, 27.7, 29.0, 28.8. This 
shows that the pattern speed of the MW bar and bulge in absolute 
units, as determined by the combination of the RCGs density and 
the BRAVA data, is between 25 and 30kms -1 kpc -1 . This is slightly 
lower than the value of £2 p ~ 30 — 40kms -1 kpc -1 determined by 
|Long et aT7| ( |2013| ), also from the BRAVA data. The comparison with 
other determinations in the literature is discussed in Section [831 
Our initial models were constructed with different dark matter 
haloes, and the bars formed in these models had different corota¬ 
tion radii and & values & = 1.1 — 1.8 based on their individual 
rotation curves (Fig.[I]). After the rescaling during the NMAGIC fit, 



the rotation curves of all models are essentially identical in the in¬ 
ner 3 kpc, such that these models all provide equally good fits to the 
RCGs density and BRAVA data. This is displayed in Fig. [9] which 
shows the azimuthally averaged rotation curves of our best mass 
models, together with their disc and halo contributions. The mod¬ 
els’ scaled outer rotation curves are different, however, consistent 
with the different values. This is not expected to influence the 
bulge dynamics, because very few stars from 4 kpc and beyond will 
reach the BRAVA bulge fields along their orbits. 

Because no effort has been made to match the rotation curve 
of the MW at large radii, the models’ ratio does not correspond 
to that of the MW. There are several ways in which the outer rota¬ 
tion curve of the MW could have changed after the bar and bulge 
formed, e.g., by later growth of the disc. Therefore, in order to es¬ 
timate the corotation radius and of the Galaxy corresponding 
to £2 p ~ 25 — 30kms -1 kpc -1 we need to use additional data. As¬ 
suming the composite rotation curve of |Sofue, Honma & Omodaka] 
(2009) rescaled to (Ro,Vq) = (8.3kpc,218kms -1 ), this range of 
pattern speed would result in a corotation radius between 7.2 and 
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Figure 9. Circular velocity curves obtained from the azimuthally averaged 
potential of our best dynamical mass models after NMAGIC fit. The solid 
lines display the total rotation curves while dashed and dotted lines display 
the halo and disc contributions, respectively. 


8.4kpc, which implies between 1.5 and 1.8. The MW would 
then belong to the so-called slow rotators. This result is discussed 
in more detail in Section IQ1 


4.4 Proper motions 

As an independent check we predict proper motion dispersions <7/^ 
in the / and b direction for our best dynamical models as explained 
in Section |3.3.3[ and compare them to the data from Rattenbury 
|et al.| ( |2007a| . A comparison model/data is shown in Fig. [lO] for 
<7/ (left-hand plot) and Oj 9 (right-hand plot). The different colours 
refer to the five best dynamical models, and the shaded regions dis¬ 
play different levels of relative error of the model with respect to 
the data. Data error bars are not plotted: errors given by [Rattenbury | 
|et al. | ( [2007a) are only statistical errors at the level of 1% or bet¬ 
ter. However, reproducing their derivation of (J/ ^ from the original 
motions of individual stars in a couple of fields, it seems to us that 
systematic errors dominate. These systematic effects are due to the 
selection thresholds and have a typical magnitude of 10%, which is 
indicated in Fig.[K)]by the white band. 

Our models provide good proper motion predictions in the l 
direction, being mostly inside the 10% limit. Proper motions in the 
b direction are slightly worse, with model values mostly 10% to 
20% lower than the data. This hints at additional systematic errors, 
either in the models or in the data, a& is directly related to the verti¬ 
cal derivative of the potential and therefore to the mass concentra¬ 
tion toward the midplane. This is apparent in Fig. [TO] where models 
with a more massive stellar bulge component have larger o^. There 
is only limited scope for increasing a& in the models as will be 
discussed in Section |0| [Rattenbury et al.| ( |2007a| ) noted that devia¬ 
tions of proper motion dispersions between adjacent fields could be 
as high as 0.2masyr -1 and concluded that some small systematic 
effect could indeed be present in the data. All together our models 
predict reasonable proper motion dispersions, indicating that the 
dynamics of the GB is quite constrained by the combination of the 
BRAVA data and the RCGs density. 
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Figure 11. Mass of the bulge-in-box for our five best mass models in dif¬ 
ferent dark matter haloes. The blue curve refers to the stellar mass while the 
red curve refers to the total mass. 


5 MASS OF THE GALACTIC BULGE 
5.1 Evaluation of the total mass 

The stellar mass in the bulge-in-box (b-b, the inner ±2.2 x ±1.4 x 
±1.2kpc of the GB) is determined by the value of & derived in 
the previous section together with the velocity scaling found by 
NMAGIC during the fit. The stellar mass recovered this way dif¬ 
fers from model to model, ranging from about 1.25 x 10 10 Mq for 
model M80 to 1.6 x 10 10 Mq for model M90. This is expected given 
the fact that our models have different dark matter masses. More 
massive haloes (like M80) can build the BRAVA dispersion with 
relatively low stellar mass while low mass haloes need more stel¬ 
lar mass. All our models are good fits to the data, so purely from 
the BRAVA data we cannot infer the stellar mass of the bulge accu¬ 
rately. However, our modelling gives us a very good estimate of the 
total mass of the bulge-in-box. Fig.[lT]shows the stellar mass (blue 
points) and total mass (red points) of the b-b for all models. We can 
see that our estimates of the total mass are quite constant along our 
range of model haloes. Altogether we evaluate the total mass of the 
b-b to be 1.84 ±0.07 x 1O 1O M 0 . The errors quoted here are com¬ 
bined statistical and systematic whose determination is discussed 
below. 

To estimate statistical errors, we use a Ax 2 analysis, based on 
two approaches. First, we regard the velocity dispersion profiles 
for a given model with different & parameters as a one-parameter 
family of curves matched to the 82 BRAVA velocity dispersions, in 
which case the appropriate A% 2 = 1. This leads to an average un¬ 
certainty for the models of A& = 0.04, AM s = 0.025 x 1O 1O M 0 , 
AM to t = 0.028 x 1O 1O M 0 . Secondly, we consider the kinematic 
Xk/ n k~v alues of the best-fitting models from Table[3]for the differ¬ 
ent dark matter haloes as the combination of a systematic variation 
modelled as linear, plus a fluctuating component which has a root 
mean square of rms(^) = 0.0119x164 =1.95. We take this as an 
approximation of the scatter in x% at minimum, which according to 
the simulations of |Morganti et al.| < [2013]) can be taken as a proxy 
for the Ax 2 to be used for estimating the accuracy with which the 
mass can be determined from the data using NMAGIC modelling. 
Applying this to the combined Xk curve derived from Fig. [^results 
in estimated uncertainties of A, & = 0.05, AM s = 0.033 x 1O 1O M 0 , 
AM to t = 0.036 x 1O 1O M 0 . Based on both methods, we estimate the 
statistical uncertainty in the stellar and total mass measurement for 
the b-b as AM s = 0.03 x 10 10 M o , AM tot = 0.03 x 1O 1O M 0 . 
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I I < 10% IZZI < 20% IZZI < 30% IZZI > 30% 



M80 I I M82.5 I I M85 I I M87.5 M90 

Figure 10. Comparison between model and data for proper motion dispersions in the l direction (left-hand plot) and in the b direction (right-hand plot) for all 
bulge fields from |Rattenbury et*aL]j2007a) . The different colours refer to the models as stated in the legend. The shaded regions indicate the relative error of 
model predictions with respect to the data, from ^ 10% in white to ^ 20%, ^ 30% and ^ 30% in dark grey. 


5.2 Evaluation of systematics 

Our modelling relies on some assumptions whose influence on the 
derived bulge mass we now investigate. Here we describe varia¬ 
tions of our four main assumptions: the midplane extrapolation, the 
length scaling, the snapshot selection and the bar angle. We show 
that the effect of these assumptions on the estimate of the mass of 
the bulge is small. These test variations are then used to set the error 
on the mass measurement previously quoted. 


5.2.1 Midplane extrapolation 

In order to quantify the uncertainty introduced by the assumed 
shape of the extrapolation in the midplane (see Section |T2) , we do 
the same study as detailed in Section[4]with a variant extrapolation, 
more concentrated in the midplane. 

In the fiducial case, we use the best sech 2 fit to fill in the 
midplane. For our variant, we perform an exponential extrapolation 
down to the midplane, with scale height h varying as a function of 
the in plane coordinates (V,/) (see Section [2^2| . We assume that 
h(x',y') is proportional to the scale height of the best sech 2 fit at 
large z, normalized with numerical value fixed such that /i(0,0) is 
equal to 1° ~ 140pc. The implied additional RCGs in the midplane 
strip increase the total number of RCGs in the b-b by 10% with 
respect to our fiducial extrapolation. We consider this extrapolation 
as giving extreme but still reasonable stellar concentration towards 
the midplane. 

The vertical density profiles along the major axis of this vari¬ 
ant extrapolation are plotted along with our fiducial ones in Fig.|T2| 
In this figure, the different colours indicate the absolute value of the 
position along the x' axis, from red at the centre to blue at the edges 
of the 3D map. 



z [kpc] z [kpc] 

Figure 12. Vertical density profiles along the major axis of the bar for the 
fiducial extrapolation of the density to the midplane (left plot) and the vari¬ 
ant extrapolation (right plot). The different colours indicate the coordinate 
along the major axis, from red in the centre to blue at the edge of the map at 
2.2kpc. The colour shaded region depicts the errors in the density from 
|Wegg & Gerhard||2013| and the dashed lines show the ±150pc region 
where the original density map from Wegg & Gerhard (2013) has been ex¬ 
trapolated. 


5.2.2 Length scaling 

In the fiducial case we scale our model using the length of the long 
bar, assuming that this long bar ends at / = 27°. Even though this 
assumption seems well founded, such a long bar is in clear ten¬ 
sion with previous claims that the pattern speed of t he M ilky Way 
bar could be as high as 60km s -1 kpc -1 (see Sectionjo]). Ongoing 
work by Wegg et al. (in preparation) shows that the long bar can 
be reliably traced up to at least / = 20°. Hence we repeat our ex¬ 
periments by scaling our models on a long bar which would end at 
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/ = 20° instead of 27°. With the assumed bar angle and distance to 
the GC the semi-major axis of such a bar would be 3.8kpc long. 


5.2.3 Snapshot selection 

Our fiducial study is based on initial models which are late evo¬ 
lutionary snapshots of pure disc+halo simulations. Throughout its 
evolution the bar gives away angular momentum, slows down and 
builds a strong B/P bulge. As our target density has a very strong 
peanut shape, late evolutionary stages are a priori more suitable 
starting points for our modelling. However, when looking at the 
ratio of the size of the peanut shape to the length of the bar, we 
found that early evolutionary snapshots better match our target ra¬ 
tio. Consequently we repeat the same modelling analysis using ear¬ 
lier snapshots, taken just after the buckling is complete. 


5.2.4 Bar angle 


The angle between the major axis of the barred bulge and the Sun- 
GC line has in the past generally been found in the range 20 — 30° 


(Stanek et al. 1997 

Bissantz & Gerhard 2002, Rattenbury et al. 

2007b, Nataf et aL 

2013). In this study we assumed an angle of 


27° which is what WG13 measured with an accuracy of ±2° when 
making their 3D density map of RCGs in the bulge. Even though 
this result seems robust we quantify the effect of a different bar 
angle by experimenting with a bar angle of 32° (2.5<r). 


5.2.5 A very robust estimate of the total mass 

Fig. [13] shows the total mass M tot as a function of the model for 
our fiducial case as well as for each variant detailed above. The 
average total mass in the fiducial case gives an estimate of the 
total mass of the bulge of 1.84 x 1O 1O M 0 , shown by the dashed 
line in Fig. Systematic variations of our four main assump¬ 
tions have only small effects on the derived value of the total 
mass as shown in Fig. [13] As stated in Section |4.1| the uncer¬ 
tainty of the total mass determination due to the discrete sam¬ 
pling of & is less than 0.01 x 1O 1O M 0 . The estimated statistical 
error is A sta xMt ot = ±0.03 x 1O 1O M 0 . Systematics dominate and 
are evaluated by simply taking the range of all mass measure¬ 
ments, as showed by the grey band in Fig. El corresponding to 
A sys M tot x 10 10 M S ~ ±0.06 x 1O 1O M 0 . Adding the statis- 

tical and systematic error in quadrature, we conclude that the total 
mass of the bulge-in-box is 1.84 ±0.07 x 1O 1O M 0 . Numerical val¬ 
ues of stellar, dark matter and total mass in the b-b for the different 
models are given in Table [4] 


6 MASS-TO-LIGHT AND MASS-TO-CLUMP RATIOS 

A good estimate of the total mass together with a constraint on the 
stellar mass would give useful insight into the dark matter mass 
in the bulge. In this section we constrain the stellar mass in the 
b-b through the stellar mass-to-light ratio T k and what we call the 
stellar “mass-to-clump” ratio, i.e. the amount of stellar mass per 
number of Red Clump and Red Giant Branch Bump stars. We first 
construct predictions from population synthesis models and com¬ 
pare our best mass dynamical models to these predictions. In this 
way we can relate the stellar Initial Mass Function (IMF) to the 
dark matter mass of the GB. 



Extrapolation 


Snapshots selection 


] Bar length 


Figure 13. Total mass (stellar + dark matter) of the bulge-in-box for the 
five models. The black points show our fiducial results already plotted in 
Fig. ED The coloured dots show masses obtained when varying some of the 
model assumptions as stated in the legend. Only small deviations from our 
fiducial case are found. The dashed line displays the mean value obtained 
in the fiducial case while the grey band span the range of all results. 



] Kroupa Zoccali 


Figure 14. Plot of the logarithmic form of the different IMF used in this 
study jSalpeter|1955[|Zoccali et al.|2000[|Kroupa|2001[|Chabrier|2003) . 


6.1 Population synthesis models 

We predict the stellar mass-to-light and mass-to-clump ratios from 
modelling the evolution of the bulge stellar population. This mod¬ 
elling relies on the three different ingredients. 

(i) An Initial Mass Function (IMF): 

We use four IMFs which span the range of reasonable IMFs for the 
Galac tic b ulge. Our two extremes are the bottom-heavy Salpeter 
IMF |Salpeter||l955| > for masses between 1O _1 M 0 and 1O 2 M 0 , 
and the Zoccali IMF ( [Zoccali et al.|2000| third entry of their table 
3), the latter being derived specifically from measurements towards 
the Galactic bulge. In between we use the Kroupa ( |Kroupa|2001) 
and the Chabrier IMF (Chabrier 2003), which are quite similar. 
These four IMFs normalized to 1M 0 are displayed in Fig.[l4] Fol¬ 
lowing [Krou^|200l]) we plotted the logarithmic form of the IMF 
^L(logio(M)) = Mln(10) where £ L (log 10 (M))dlogi 0 (M) cor- 
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Table 4. Stellar mass M s , dark matter mass Mdm and total mass M tot for all models in the bulge-in-box, in units of 1O 1O M 0 . The first row refers to our fiducial 
models while the lower four rows refer to the different systematic variations of the model assumptions detailed in Section [52] 



M s 

M80 

Mdm 


M s 

M82.5 

Mdm 


M s 

M85 

Mdm 


M s 

M87.5 

Mdm 


M s 

M90 

Mdm 


Mtot 

M to t 

M to t 

M to t 

M to t 

Fiducial 

1.27 

0.53 

1.80 

1.37 

0.45 

1.82 

1.47 

0.38 

1.85 

1.59 

0.29 

1.88 

1.63 

0.23 

1.85 

Extrapolation 

1.29 

0.54 

1.83 

1.36 

0.47 

1.83 

1.49 

0.37 

1.86 

1.58 

0.29 

1.87 

1.63 

0.23 

1.86 

Bar length 

1.32 

0.55 

1.87 

1.41 

0.50 

1.91 

1.50 

0.38 

1.88 

1.57 

0.30 

1.87 

1.60 

0.24 

1.84 

Bar angle 

1.30 

0.55 

1.85 

1.38 

0.48 

1.86 

1.50 

0.38 

1.88 

1.60 

0.29 

1.89 

1.64 

0.23 

1.87 

Snapshots selection 

1.41 

0.42 

1.83 

1.48 

0.37 

1.86 

1.56 

0.30 

1.86 

1.69 

0.22 

1.91 

1.70 

0.17 

1.88 


responds to the fraction of stars with mass between logio(M) and 
logio (M) +dlogio(M) and M is expressed in M 0 . 

(ii) A set of isochrones: 

We choose the solar metallicity and a-enhanced BaSTI isochrones 
( [Pietrinferni et al.||2004| ) and assume a single age population of 
lOGyr. As shown later the choice of the age has a small effect and 
is therefore not critical. 

(iii) A way to treat stellar remnants: 

Stars that evolve beyond their isochrones have to be turned properly 
into white dwarfs, neutron stars or black holes. We use the choices 
described in |Maraston| ( |1998| ). 

Using these three ingredients, we first construct a luminosity 
function <h(M^) in units of mag -1 Mq 1 by evolving 1M 0 through 
the set of isochrones for a given age, according to the consid¬ 
ered IMF mass distribution. After renormalization of the luminosity 
function to a remaining mass of 1M 0 (evolved stellar population + 
stellar remnants), the synthesized stellar mass-to-light ratio in the 
K-band T k is given by 

T k=(^J ®{Mk) io-oMMk-Mkq ) dMic j (20) 

We checked our mass-to-light ratio predictions by reproducing the 
work of |Maraston| ( |1998| > and |Percival et al.| ( |2009| > and found a 
very good agreement with both of them for the old population con¬ 
sidered here. 

In order to compute the stellar mass-to-clump ratio we use 
the same technique as in WG13. We fit an exponential background 
plus two Gaussians to the previously derived luminosity function. 
The two Gaussians represent the RCGs and the Red Giant Branch 
Bump (RGBB), as in WG13. This mass-to-clump ratio includes 
RGBB stars as well for two reasons. First because for some ages, 
the RCGs and RGBB have the same luminosity and are therefore 
indistinguishable. Secondly because WG13 made their map by fit¬ 
ting these two Gaussians under the assumption that the RGBB rep¬ 
resents 20% of the Red Clump. Hence the total number of RCGs 
and RGBB stars is a priori better constrained than the number of 
RCGs only. As the luminosity function was renormalized to a re¬ 
maining mass of 1M 0 , the number of stars contained in the two 
fitted Gaussians directly leads to the mass-to-clump ratio, denoted 
as M/n RC+ RGBB- 

6.2 Mass-to-light ratio 

The COBE/DIRBE instrument provides us with ^f-band measure¬ 
ments in many bulges fields. Here we use the data from jDrimmel &| 
Spergel ( 2001) who removed point sources by applying a median 
filter to the original “Zodi-Subtracted Mission Average” map |Kel-| 
|sall et al.[1998| ). In order to correct for Galactic extinction, we use 


the extinction map presented in WG13, also derived from K-band 
data. This map covers the inner / G [—10°,5°],Z? G [—10°, 10°] of 
the bulge and has a resolution of l' which is much finer than the 
42' x 42' COBE/DIRBE field size. We correct the COBE/DIRBE 
data by using the mean extinction on each COBE/DIRBE pointing. 
We ignore the fields in the inner \b\ <2° because the extinction is 
too high to make a reliable correction. We apply a disc contamina¬ 
tion correction to the remaining data points as follows. The average 
disc contamination as a function of latitude is evaluated using the 
surface brightness profile in two 0.5° wide strips along the b di¬ 
rection, located at / = ±15°. At these large longitudes the Galactic 
bulge is no longer important and the flat long bar does not con¬ 
tribute significantly to the surface brightness for \b\ ^ 2°. Hence 
the average surface brightness profile of the strips is mostly due 
to the disc component. By assuming that this disc vertical surface 
brightness profile stays constant at all longitudes inside |/| < 15° 
we can estimate and then remove the disc contamination from the 
COBE data. Finally this provides us with about 2800 extinction and 
foreground corrected surface brightness measurements towards the 
Galactic bulge which can be compared to our models. 

To compute the surface brightness of our particle models we 
convert scaled stellar particle weights w* into ^-band luminosity 
Li = WiYg 1 , assuming a constant but still unknown stellar mass to- 
light ratio T k in the bulge. will be determined by matching our 
model surface brightness to the COBE/DIRBE K -band data. The 
extinction-free model surface brightness in some field Zy (/,/?) as 
one would see from the Sun’s location is given by 

(21) 

J i i 

where A Llj is the solid angle of the considered field, r ; is the dis¬ 
tance of the particle i from the Sun in kpc, and Sjfc) is a suit¬ 
able selection function that we describe below in Eq.[22] Our mod¬ 
els have been matched to the MW inside the b-b but no attempts 
have been made to match the model discs to the MW disc. In order 
to reduce the uncertainty due to the disc contribution we evalu¬ 
ate for each COBE field the fraction of model light which comes 
from particles located in the b-b. We remove from the analysis all 
fields where this fraction is lower than 90%, i.e. all fields where 
not directly constrained particles contribute more than 10% of the 
model surface brightness. Given that a single particle can theoret¬ 
ically dominate the surface brightness by being arbitrarily close to 
the Sun’s location, we also remove nearby disc particles from the 
analysis in the remaining fields by taking the following selection 
function: 

1 if i G field / and n ^ 3 kpc 

F (22) 

0 otherwise 
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I I Salpeter I I Chabrier 
I I Kroupa Zoccali 

Figure 15. Stellar mass-to-light ratio in the K band for our five best dynam¬ 
ical models of the bulge in different dark matter haloes. The model errors 
shown are dominated by systematic effects. The different coloured lines 
indicate predictions from different IMFs as stated in the legend. The black 
dashed line is an estimate of the highest allowed mass-to-light ratio obtained 
by turning all dark matter of the b-b into luminous matter. This allows us to 
rule out the Salpeter IMF for the Galactic bulge with age lOGyr. 


We then apply the same foreground contamination correction to 
the model surface brightness as for the data, using the two strips 
at / = ±15°. We checked that our results do not depend on the 
exact form of the selection function, indicating that the foreground 
contamination has been properly removed. 

Finally, we compute the mass-to-light ratio as stated in Eq.[2T] 
independently for all remaining COBE fields and average the re¬ 
sults. The statistical error in the mean mass-to-light ratio is very 
low due to the large number of COBE fields so that systematics 
dominate. We evaluate systematic effects by repeating this analy¬ 
sis for the four variants described in Section [5^21 of the considered 
model. The full range of values is then taken as the systematic error. 

Fig. |T5] shows the mean value of T k and its associated sys¬ 
tematic error for the five best dynamical models with different dark 
matter haloes. These mass-to-light ratios are compared to the pre¬ 
dictions from different IMFs shown as the coloured lines, for a sin¬ 
gle lOGyr age population. The coloured strips span the range of 
values for ages between 9 and 11 Gyr. With the total mass and inte¬ 
grated light fixed, the models with low dark matter mass (M90 for 
example) have higher stellar mass-to-light ratios. We see that most 
of our models are in agreement with predictions from a Kroupa or 
Chabrier IMF which are somewhat similar, the Kroupa IMF pre¬ 
dicting slightly higher mass-to-light ratios than the Chabrier IMF. 
Only the most dark matter dominated model M80 approximately 
matches the predictions from the Zoccali IMF, which is a priori 
the best candidate IMF because it was measured directly from the 
stellar luminosity function of the bulge. In order to agree with the 
Zoccali IMF about 40% dark matter mass is required in the b-b. 

In addition, the dynamical models rule out a Salpeter IMF 
for a bulge population with age lOGyr. The dashed line in Fig. 
p3] represent the highest possible stellar mass-to-light ratio which 
would be obtained if all the dark matter in the bulge was turned 
into stars with the same density distribution as the stellar compo¬ 
nent of the bulge in our fiducial models. For all our models, the 
mass-to-light predictions for a Salpeter IMF for a lOGyr Galactic 
Bulge are higher than this extreme mass-to-light ratio by at least 
three times the model error, showing that it is too bottom-heavy. 



Figure 16. Mass-to-clump ratio, M/nRc+RGBB for our five best dynamical 
models of the bulge in different dark matter haloes. The coloured lines in¬ 
dicate predictions from different IMFs as stated in the legend. The black 
dashed line is an estimate of the highest allowed mass-to-clump ratio ob¬ 
tained by turning all dark matter of the b-b into luminous matter. 


6.3 Mass-to-clump ratio 

The mass-to-clump ratio is also a useful quantity to relate stellar 
mass and stellar population. With RCGs being approximate stan¬ 
dard candles with standard colours, issues like foreground con¬ 
tamination and dust extinction are easier to solve when computing 
a mass-to-clump ratio than a mass-to-light ratio. The number of 
RCGs in the b-b was computed directly from our fiducial extrap¬ 
olation of the density map of WG13. Using the variant extrapola¬ 
tion presented in Section |5.2| would increase this number by only 
10%. In order to include the RGBB in this calculation, we add the 
20% fraction assumed in the derivation of the 3D map by WG13. 
The computation of M/nRc+RGBB is then straightforward from the 
stellar mass determination of Section[5J] Results are shown in Fig. 
p~6|for all dark matter models. The errors plotted in Fig.[l6]are sys¬ 
tematic, determined from the four variant assumptions described in 
Section |5^2| Again the different colour strips indicate predictions of 
different IMFs for ages between 9 and 11 Gyr. 

We reach similar conclusions from the mass-to-clump ratio as 
from the mass-to-light ratio. This indicates that the issues of ex¬ 
tinction and foreground contamination present in the computation 
of the mass-to-light ratio have been treated correctly. A 35% dark 
matter fraction in the b-b is needed to match predictions from the 
Zoccali IMF and low dark matter models would need a Chabrier 
IMF to match the predictions for a lOGyr old population. Once 
again the Salpeter IMF over-predicts the mass-to-clump ratio by 
about 50% and can therefore be ruled out for a lOGyr GB. 


7 PEANUT AND X-SHAPE STRUCTURES OF THE 
GALACTIC BULGE 

In this section we study the structural properties of the GB further, 
based on the 3D density map of RCGs from WG13. We first il¬ 
lustrate and discuss its X-shape in a similar way as usually done 
for external galaxies. We then use a 3D photometric diagnostic to 
quantify the fraction of stellar mass involved in the peanut shape of 
the GB. 
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Figure 17. X-shape structure of the Galactic Bulge obtained by unsharp 
masking as described in |Bureau et al.|p006) . The Galactic bulge has a so- 
called off-centred X-shape. 

7.1 The photometric X-shape 

Observations of external galaxies have revealed that Box/Peanut 
bulges exist with a variety of shapes. In order to highlight the in¬ 
ternal structure of external edge-on B/P bulges, a common practice 
is to use unsharp masking techniques as described by Bureau et al. 
P006| . They applied a median filter to images of 30 edge-on spi¬ 
rals and removed the smoothed images from the original ones. This 
reveals what is called the X-shape. Bureau et al.|( |2006| > proposed a 
classification of external B/P bulges based on the properties of this 
X-shape. Where would the Galactic bulge appear in such a clas¬ 
sification? To answer this question, we apply a median filter with 
kernel size 500pc to the side-on projection of the 3D density of 
RCGs from WG13 and remove it from the original projection. The 
positive residuals revealing the X-shape of the GB are shown in 
Fig. [TT] with contours spaced by a third of a magnitude. 

Fig. [TT] shows that the Galactic bulge has an off-centred X- 
shape structure, its two arms crossing the major axis about 500pc 
away from the centre. This feature is identified in 50% of external 
edge-on B/P bulges ( [Bureau et al.|2006| ). 

7.2 The mass of the peanut shape 

Consistent with the strong X-shape, the Milky Way bulge in the 
side-on map from WG13 shows a very prominent peanut shape. An 
interesting issue is the amount of stellar mass involved in this fea¬ 
ture. This was already addressed by [Li & Shen|| ( [2012| who applied 
a technique somewhat similar to unsharp masking to the side-on 
projection of the model of |Shen et al.| ( |2010] ). They fitted ellipses to 
the isophotes of the side-on projection, modelled the light projected 
by these elliptical isophotes and removed the modelled light from 
the original image. Doing so revealed a centred X-structure in the 
model accounting for about 7% of the bulge stellar mass. Because 
their model had been shown to give a good representation of the 
BRAVA kinematic data, they then concluded that the stellar mass 
involved in the peanut shape of the GB was similarly about 7%. 

Since we now have a direct measurement of the 3D density 
of RCGs in the bulge, we can perform a similar photometric anal¬ 
ysis directly on the 3D density map. Consider the density profile 
p (0,0, z ) along the minor axis of the bar. For each particular value 
of z, we evaluate the 3D isodensity surface with density p(0,0,z). 
Then we look for the most voluminous ellipsoid we can find which 
fits inside this particular 3D isodensity surface. This search is done 
under the constraint that the ellipsoid is centred on the centre of the 


£j?c[pc] 2 

0.10 0.24 0.56 1.33 3.16 




Figure 18. “Non-ellipsoidal” residual density of RCGs in the density map 
of |Wegg & Gerhard] (2013) . The residual density accounts for 24% of the 
RCGs in the bulge. Plotting conventions are the same as in Fig. [3] 


bulge, the principal axes are aligned with the principal axes of the 
bulge, the semi-principal length along the vertical axis is fixed to 
|z|, and the ellipsoid stays inside the isodensity surface considered. 
By doing so for all z we construct a family of ellipsoidal isoden¬ 
sity surfaces. We compute the 3D density arising from this ellip¬ 
soidal isodensity family and remove it from the original 3D map. 
The residuals correspond to the “non-ellipsoidal” part of the bulge 
density. As all our ellipsoids are truly inside their corresponding 
isodensity surfaces of the original map, we are assured that the 
residual map is positive at all points. 

To reduce the effect of measurement errors in the original map 
of WG13, we actually do this analysis using model M80, which 
through NMAGIC fitting of the density gives a smoother density 
map which is everywhere within the errors of the original map. We 
note that the results do not depend on which initial model is used 
for the density fit. 

The residual map is plotted in projection in Fig.[l8] The four 
lobes responsible for the peanut shape of the bulge are clearly visi¬ 
ble in the side-on view for \z\ > 300pc. For w < 300 pc the peanut 
is not prominent enough with respect to its surroundings to make 
the density shape deviate from ellipsoidal shape. The stellar mass 
involved in this residual peanut shape is about 24% of the total 
stellar mass of the bulge. This figure probably underestimates the 
real amount of mass involved in the peanut structure as one expects 
the orbits responsible for this structure to also visit the strip inside 
\z\ < 300pc and therefore to contribute to the ellipsoidal density 
that was removed here. 

We can estimate an error on the 24% by applying the same 
procedure to each of the modelled five variant maps of WG13. We 
find that the non-ellipsoidal residuals account for 24^% of the 
stellar mass of the bulge. By applying our diagnostic in a two di¬ 
mensional version to the side-on projection of the density, we find 
that the 2D residuals sum to only 1l+}% of the stellar mass of the 
bulge, showing that 2D determinations tend to underestimate the 
mass in the peanut structure. 

Our peanut mass fraction of 24^% results from an observa¬ 
tional definition of the peanut shape: a deviation from ellipsoidal 
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density shape. A more physical definition of the peanut structure 
would be to identify the different orbits responsible for its shape 
in side-on projection. We are currently working on such an orbit- 
based characterization of the peanut shape which will be presented 
in a companion paper. 


8 DISCUSSION 


8.1 Dynamical mass of the bulge 

The Galactic B/P bulge transits into a longer two-dimensional bar; 
therefore in this paper we use a simple definition for determining 
the dynamical mass of the bulge. We use the three-dimensional box 
±2.2 x ±1.4 x ±1.2kpc as our bulge volume in which the RCGs 
density was determined by Wegg & Gerhard (2013) and which con¬ 
tains most of the three-dimensional bulge part of the bar. The total 
mass of this bulge-in-box (b-b) is accurately determined by the dy¬ 
namical models, M tot = 1.84 ±0.07 x 1O 1O M 0 . This value is es¬ 
sentially independent of the dark matter mass fraction in the bulge 
throughout our models, and the estimated error includes both sta¬ 
tistical and systematic uncertainties ; the systematic part is deter¬ 
mined from varying the modelling assumptions around our fiducial 
model. 

Dynamical models have previously been used to estimate the 
mass of the bulge. Using the Schwarzschild method, [Zhao| ( | 1994) 
built a self-consistent model of the bar/bulge and found a total 


bulge mass of 2 x 1O 1O M 0 while Kent (1992 ) found a mass of 


1.8 x 1O 1O M 0 by modelling an oblate isotropic rotator with a con¬ 
stant mass-to-light ratio. By studying gas dynamics in the poten¬ 
tial of the model from Bissantz & Gerhard ( 2 002|> obtained by de- 
projecting the COBE luminosity distribution, |Bissantz, Englmaier| 
|& Ge rhard (2003) determined the circular velocity at 2.2kpc to 
be 190kms . Once converted to mass under the assumption of 
spherical symmetry, this leads to a total bulge mass of about 
1.85 x 1O 1O M 0 . All these results compare well with our estimate 
of the mass of 1.84 ±0.07 x 1O 1O M 0 , especially when considering 
the difficulty of a precise definition of the bulge. 

An independent way to obtain a dynamical mass is from the 
virial theorem. Such studies lead to quite different values depend¬ 
ing on the assumed p attern speed and bar angl e, from 1.6 x 1O 1O M 0 
jHan & Gould|l995) , up to 2.8 x 1O 1O M 0 jBlum|l995| obtained 
for a high pattern speed (81 kms -1 kpc -1 ). The use of the virial 
theorem has two weaknesses. First, it relies on an estimation of 
the total kinetic energy, which cannot be measured accurately from 
line-of-sight data only. Secondly, it assumes that the bulge on its 
own is a system in equilibrium, whereas it is in fact part of a bar 
embedded in a disc. 

Traditionally, the contribution of dark matter to the mass in 
the bulge region has been considered unimportant. In this case, a 
bulge mass can be estimated simply from the stars. For example, 
|Dwek et ak] ( |1995| estimated the stellar mass of the bulge at 1.3 x 
1O 1U M 0 from the COBE luminosity and a Salpeter IMF. However, 
neglecting dark matter in the bulge is not a good approximation as 
we have seen in Section [5] 


8.2 Stellar and dark matter mass in the bulge 

Because the total mass of the bulge-in-box (b-b) is well determined, 
knowing the stellar mass independently would give insight into the 
amount of dark matter in the central part of the Galaxy. As shown in 


Section [6] the stellar mass can be constrained by comparing mass- 
to-light and mass-to-clump ratio measurements with stellar pop¬ 
ulation synthesis predictions. However, for this the choice of the 
IMF is crucial: candidate IMFs disagree within a factor of two in 
their prediction of the stellar mass-to-light ratio, so a reliable mea¬ 
surement of the IMF in the GB is needed. Currently, the Zoccali 
IMF ( [Zoccali et al.|2000) measured from the bulge star luminosity 
function in a field located near (l,b) = (0°,—6°) is the favoured 
IMF for the GB. For this case, the predicted mass-to-light ratio 
from isochrones (Section |6.1| ) gives a stellar mass for the b-b of 
1.12 x 1O 1O M 0 . From our modelling, both the mass-to-light and the 
mass-to-clump ratio independently agree that models with fairly 
high dark matter mass fraction are required to provide the remain¬ 
ing part of the dynamical mass in the b-b. These models predict a 
dark matter mass in the b-b of ^ 0.7 x 1O 1o M 0 , accounting 
for about 40% of the total mass of the b-b. 

Further insight on the dark matter part of the b-b mass can be 
obtained from proper motions. In particular, the proper motion dis¬ 
persions in the b direction directly constrain the derivative of the 
potential along the vertical direction, and therefore the total mass 
concentration towards the plane. We showed in Section |4.4| that 
our proper motion predictions are 10% to 20% lower than the 
data. Increasing the stellar mass concentration towards the plane as 
in Section [5.2.1| can indeed increase o^, but only at percent level 
which is not significant enough. It is also possible that a systematic 
effect is present in the data, e.g., due to extinction or incomplete¬ 
ness of the sample. Before we can use the proper motion data to 
measure the mass concentration and dark matter content of the GB, 
these possible systematic effects in the data need to be better un¬ 
derstood. A more detailed study of the proper motion constraints is 
part of our ongoing work to make a more complete model of the 
Galactic bulge and long bar. 


8.3 Pattern speed of the MW bar/bulge 


Our dynamical models, based on the RCGs density from |Wegg| 
|& Gerhard] ( [2013) ) together with the BRAVA kinematic data for 
the bulge stars, imply a pattern s peed of the MW bar-bulge of 
25 — 30kms -1 kpc -1 (see Section [±3| . This result remains un¬ 
changed when varying the modelling assumptions, for example 
when we consider a shorter bar, as detailed in Section |5.2| Com¬ 
paring with the composite rotation curve from |Sofue, Honma &| 
|Omod aka ( 2009]), this would place the corotation value of the bar 
just inside the solar circle, and give a ratio of corotation over bar 
half length between 1.5 and 1.8. The MW would then belong to 
the so-called slow rotators. Slow rotators are quite rare for external 
galaxies in general ( |Aguerri, Beckman & Prieto|1998|), but fairly 
common among SBc galaxies (Rautiainen, Salo & Laurikainen| 
j 2008 |. 

There have been a quite a number of pattern speed measure¬ 
ments in the MW by other techniques. The bulk of the measure¬ 
ments indicate a fast bar, but some studies suggest a lower pat¬ 
tern speed (se e|Gerhard|p011|) fo r a review). Direct measurement 
by |Debattista, G erhard & Sevenster (20 02]) using a variant of the 
Tremaine & Weinberg method ( [Tremaine & Weinberg) 1984f leads 
to the high value of O p = 59 ±5kms -1 kpc -1 but depends strongly 
on the radial velocity of the local standard of rest towards the GC. 

Numerous indirect measurements of the bar pattern speed 
have been obtained by matching some of the observed fea¬ 
tures in the position-velocity diagrams for Hi and CO to gas- 
dynamical model predictions. Most of these studies argue for 
a high pattern speed, such as |Englmaier & Gerhard| \\999) 
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who found ~ 60kms ^pc 


.-l 


Fux 


jl999t who obtaine d 


50kms 1 kpc 1 , and Bissantz, Englmaier & Gerhard (2003) with 


55 —65kms ^pc ^However, 


Weiner & Sell wood ( |1999| > and 
( |2008| ) obtained lower values of 


Rodriguez-Fernandez & Combes 
respectively 42km s -1 kpcand 30 — 40kms -1 kpc -1 . 

High pattern speeds were also obtained by explaining the stel¬ 
lar kinematics of nearby disc stars with the dynamical effects of the 
Outer Lindblad Resonance (OLR) of the bar. By interpreting the bi¬ 
modality of th e velocity distribution in the solar neighbourhood in 
this way, 


Dehnen 


2000 found O p = 53 ± 3 kms 1 kpc 1 . Antoja 


|et al. | ( |2014[ ) presented an analytical model of the effect of the OLR 
in the velocity distributions of stars at different radii and showed 
that they could reproduce measurements of the Hercules stream for 
a bar pattern speed of O p = 54 ± 0.5 kms -1 kpc -1 . 

The low pattern speed value found from the RCGs and BRAVA 
data is consistent with |Rodriguez-Fernandez & Combes| ( [2008) and 
|Long et ak] ( |2013| but not with the majority of these measurements. 
We note that models based on variants of the RCGs density map 
as well as models constrained by a model B/P bulge density inde¬ 
pendent of the RCGs measurement jMartinez-Valpuesta[2012| give 
similarly low values. The BRAVA kinematics are in agreement with 
the ARGOS kinematics ( |Ness et al.|2013| despite the quite different 
selection functions of both surveys. In our modelling, the details 
of the selection function assumed for the BRAVA data were not im¬ 
portant. Therefore our value measured from modelling the bulge 
kinematics appears robust. 

The highest pattern speeds reported in the literature are also 
in conflict with other data. Using the composite rotation curve 
from|Sofue, Honma & Omodaka|(|2009|, a pattern speed of £2 p = 
55kms -1 kpc -i would place corotation at 3.7kpc. Because the 
long bar cannot extend beyond corotation jContopoulos|1980| , this 
is incompatible with the apparent end of the long bar at / = 27° 
( |Hammersley et al.|2000[|Cabrera-Lavers et al.|2008| ). Even if we 
relax the hypothesis that the long bar ends at / = 27°, recent star 
counts from the UKIDSS survey (Wegg et al. in preparation) show 
that it can be reliably traceable at least out to / = 20°. This gives a 
lower bound of the half-length of the bar of 3.8kpc, still in tension 
with a pattern speed of 55 km s -1 kpc -1 . 

We conclude that despite of the many effort made by different 
groups, the question of the pattern speed of the Galactic bar remains 
an unsolved issue. One way to settle this question by dynamical 
modelling would be to include data that more accurately constrain 
the long bar component, which is one of our future goals. 


9 CONCLUSION 

In this work we have presented a set of self-consistent dynami¬ 
cal models of the Galactic bulge with different dark matter haloes, 
which match recent data on the spatial distribution and kinematics 
of bulge stars. We started with a family of N-body models of barred 
discs with B/P bulges, evolved from near-equilibrium stellar discs 
embedded in live dark matter haloes. We then fitted these models 
to the recent 3D density measurements of Red Clump Giants in the 
bulge from |Wegg & Gerhard) (|2013} , as well as to the BRAVA kine¬ 
matic data of Kund er et al. | ( |20 1 2] ) in multiple bulge fields, using 
the NMAGIC Made-to-Measure method. 

From this modelling, we obtain an accurate and robust esti¬ 
mate of the total mass (stellar and dark matter) of the bulge in the 
RCGs box, of M tot = 1.84 ±0.07 x 1O 1O M 0 . We also find a low 
pattern speed of about 25 — 30km s -1 kpc -1 , which with the mea- 
sured rotation curve places the Milky Way among the slow rota¬ 


tors {& ^ 1.5). We compute the mass-to-light and mass-to clump 
ratios and compare them with theoretical predictions from popu¬ 
lation synthesis models using different IMFs. We show that the 
Salpeter IMF jSalpeter[|1955| is ruled out for a lOGyr old bulge 
population. We find that a relatively high dark matter mass frac¬ 
tion in the bulge is needed in order to match predictions from the 
IMF inferred from the stellar luminosity function in the upper bulge 
|Zoccali et al. 2000), ~ 40% or Mdm ~ 0.7 x 1O 1O M 0 . In addition 
we study the X-shape of the Galactic bulge and find an off-centred 
X-shape, common in external B/P bulges ( [Bureau et al.|2006) . By a 
three-dimensional analysis of the isodensity surfaces of the RCGs 
density we find that the peanut-shaped deviations from ellipsoidal 
shape account for 24^% of the bulge stellar mass, significantly 
larger than the previous estimate of 7% of |Li & Shen| ( |20T2 ). 
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